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Equation-Free Particle-Based Computations: Coarse Projective 
Integration and Coarse Dynamic Renormalization in 2D 
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Equation-free approaches have been proposed in recent years for the computational study of 
multiscale phenomena in engineering problems where evolution equations for the coarse-grained, 
system-level behavior are not explicitly available. In this paper we study the dynamics of a diffusive 
particle system in a laminar shear flow, described by a two-dimensional Brownian motion; in partic- 
ular, we perform coarse projective integration and demonstrate the particle-based computation of 
coarse self-similar and asymptotically self-similar solutions for this problem. We use marginal and 
conditional Inverse Cumulative Distribution Functions (ICDFs) as the macroscopic observables of 
the evolving particle distribution. 
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1. Introduction 

Multiscale phenomena arise naturally in science and 
engineering. The ability to properly resolve such phe- 
nomena and propagate their influence across scales un- 
derpins the predictive value of mathematical and physics- 
based models. In the case of multiscale systems for whose 
macroscopic behavior no explicit coarse-grained, macro- 
scopic equations are available, a computer-assisted ap- 
proach, referred to as the Equation-Free Framewortii^ 
has been recently proposed. Equation-free methods nu- 
merically evolve the coarse-scale behavior through appro- 
priately designed short computational experiments per- 
formed by the fine-scale (microscopic, stochastic, agent- 
based) models. In this paper we will demonstrate the 
use of two such methods: Coarse Projective Integration 
(CPI) and Coarse Dynamic Renormalization (CDR). 

Particle-based simulators are the fine-scale description 
of choice for a variety of problems exhibiting multiscale 
behavior; such problems range from Stokesian and Brow- 
nian dynamics to the Monte Carlo modeling of microor- 
ganism locomotion, the mixing of passive scalars by tur- 
bulent velocity fields and even particle filtering applica- 
tions. The purpose of Coarse Projective Integration is to 
accelerate the computational evolution of coarse-grained 
observables of microscopic simulators; it has been suc- 
cessfully used in the past to accelerate computations of 
the collective evolution of spatially one- dimensional ran- 
dom particle distributions^. In those examples, the first 
motivated by hydrodynamics and the second by bacterial 
chemotaxis, the coarse-scale observable was the cumu- 
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lative distribution function (CDF) of the particle posi- 
tions, synthesized from snapshots of the fine-scale simu- 
lation. For one-dimensional problems in space, the func- 
tional inverse of the CDF (ICDF) is projected onto a 
suitable basis set consisting of orthonormal polynomials^ 
or POD (Proper Orthogonal Decomposition) modes (ob- 
tained through Singular Value Decomposition (SVD) of 
the sorted particle positional). The particle-level model 
is used to evolve the ICDF, and short time series of the 
coefficients of its projection on the appropriate basis are 
thus collected; these short time series are used to es- 
timate the time derivatives of the coefficient evolution. 
These estimates are then used in the context of tradi- 
tional continuum integration algorithms (such as Euler— , 
or Adams-BashforthSi) to "project" the coefficients for- 
ward in time, into the future (thus the term "projective 
integration"). To repeat the procedure, a "lifting" step 
is required: fine-scale states, i.e., particle positions whose 
ICDF is consistent with the projected coefficient values, 
are generated. Since detailed microscopic evolution over 
the duration of the projective step has been avoided, the 
procedure has the potential to alleviate the burden of full 
fine-scale simulation. Preliminary discussions on the sta- 
bility and accuracy of these schemes can be found in^. 

Beyond CPI, equation-free computational protocols 
can also be used in coarse-grained fixed point and bi- 
furcation analyses, to compute stationary states of the 
coarse-grained system dynamics and their parametric de- 
pendence. In these analyses, the action of operators on 
the coarse-grained observables is deduced from appropri- 
ately initialized computational experiments with the fine 
scale models. Matrix-free implementations of contrac- 
tion mappings, like Newton's method, have been used to 
compute fixed points of unavailable coarse-scale models 
for kinetic Monte Carlo {c.g.^, epidemiology), Brown- 
ian dynamics^ or Lattice Boltzmann^iiS fine-scale sim- 
ulators. The same approach can be used to evolve ef- 
fective medium (homogenized) descriptions of reaction- 
transport problemsii*^ based on short bursts of finely 
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For multiscale systems exhibiting scale invariance at 
the macroscopic level, renormalization techniques can 
be used to solve for self-similar solutions and their 
scalingsi^, see alsoi^*i^. Recently, dynamic renor- 
malization (e.g.^SiiiiS) was used in conjuction with 
equation-free computation to obtain coarse-grained self- 
similar solutions using short bursts of fine-level, di- 
rect simulationpiS. This coarse dynamic renormalization 
(CDR) method finds macroscopically self-similar solu- 
tions with the help of template junction^^iSli^^i^. Ob- 
serving the macroscopic solutions in a co-expanding (or 
co-collapsing) frame of reference, we seek steady states in 
the new frame; fixed point equation-free algorithms can 
be used for this task. 

Studying the evolution of particle distributions us- 
ing their ICDF as an observable is convenient for one- 
dimensional problems in space, where suitable bases for 
representing one-dimensional monotonic curves over fi- 
nite one-dimensional domains are readily available (see 
the examples in2*^iiS). The extension to corresponding 
obscrvables in more than one dimension, however, is non- 
trivial: the CDF itself, not being a bijcctivc mapping, 
does not have an inverse. In this case, operations on 
the CDF may be implemented by identifying a suitable 
set of basis functions for two or three-dimensional CDFs. 
This requires finding multidimensional orthonormal poly- 
nomial approximations for monotonic bounded functions 
with infinite support. In this paper we use an alternative 
approach, representing a multidimensional CDF in terms 
of its marginal and one-dimensional conditional distribu- 
tions. In this manner, multidimensional problems are 
converted into a collection of one-dimensional problems 
the solution of which can be obtained using standard 
approaches. Equation- free algorithms such as CPI and 
CDR can thus readily be extended to problems involv- 
ing multidimensional coarse-grained observables. Prelim- 
inary results for the self-similar case have been reported 
elsewhere24. 

The paper is organized as follows: In Section 2 a coarse 
time-stepper is constructed in terms of the marginal and 
conditional ICDFs in multidimensional particle systems. 
Use of this time-stepper in CPI and CDR is formulated in 
Sections 3 and 4, respectively. Our illustrative Brownian 
particle system in a Couette flow is described in Sec- 
tion 5, and its analytical self-similar and asymptotically 
self-similar solutions are presented. These two cases are 
then used in Section 6 to illustrate equation-free compu- 
tational procedures, and the direct, particle-based com- 
putational results are compared with the analytical solu- 
tions. We conclude with a brief summary and discussion 
in Section 7. 

2. A Coarse Time-Stepper for Multidimensional 
Random Particle Systems 

Short computational experiments with the fine-scale 
model are used to construct the coarse time-stepper ~ 
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FIG. 1: A schematic of the coarse time-stepper. 

between coarse-scale model states and fine-scale states. A 
coarse time-stepper consists of three components: lifting, 
fine-scale evolution and restriction^ (Fig- The lifting 
transformation converts coarse-scale observables to con- 
sistent fine-scale states; restriction is the reverse trans- 
formation, from fine-scale sates to coarse-grained observ- 
ables. Lifting, followed by restriction, should then give 
the identity on the coarse observables (modulo round- 
off error). Different coarse time-steppers are generated 
via different lifting and restriction operators; one should 
test on-line that macroscopic computational results are 
insensitive to the specific details of the time-stepper im- 
plementation; a more extensive discussion can be found 
ini. 

For many multiscale problems observed along a sin- 
gle effective spatial dimension, particle positions consti- 
tute the fine-scale model state, while an obvious coarse- 
grained state is the local mean concentration of the 
particles^. For identical particles this local mean con- 
centration is observed in terms of the histogram of 
the single particle position probability density function 
(PDF). However, this PDF histogram depends on the 
bin size used to estimate it; the PDF at any given point 
becomes zero if the bin size is too small so that no parti- 
cles exist within the bin containing this point. To over- 
come this difficulty, the cumulative distribution function 
(CDF) is naturally used as an alternative coarse-scale 
observable. The CDF has in principle infinite support, 
but its inverse, the Inverse CDF (ICDF), is supported in 
[0,1]; it can be readily represented by its projection on 
shifted Legendre polynomials^^. 

In multidimensional random particle systems, the in- 
verse form of the multidimensional CDF is not read- 
ily available (strictly speaking, it is not even defined). 
The marginal and conditional ICDFs constitute candi- 
date coarse-scale observables for multidimensional sys- 
tems; this is because the multidimensioiral CDF can be 
represented using the marginal CDF in one direction and 
a collection of conditional CDFs in the remaining direc- 
tions. For instance, for systems in two spatial dimen- 
sions. 



FxY{x,y) 



dF 

Fx\Y{Av)^{y)dy, (1) 



where FxYix^y), Fyiy) and Fx\Y{x\y) arc the CDF, dif- 
ferentiable marginal CDF and conditional CDFs, respec- 
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Fx\Y[x\y) = lim P{X <x\y <Y <y + /:^y) 

= lim ^xy(x-,y + Ay)-F^H.,y) 
Ay^o Friy + Ay) - Friy) 

Assuming smoothness, a finite number of conditional 
CDFs can be used to recover the particle distribution 
(e.g. through interpolation). In the following, we illus- 
trate through a two-dimensional system (without loss of 
generality) our implementation of this procedure for mul- 
tidimensional systems. 

2.1. Lifting. Starting with the inverse CDFs (ICDF) 
for the marginal and conditional distributions, the lifting 
procedure involves obtaining compatible realizations of 
the fine-scale states. Let the marginal ICDF in direction 
y, IFy{-) : [0, 1] > IR, and conditional ICDFs in the other 
direction x, IFx\y{', u) '■ [0, 1] i-^ M, y £ K, be defined by 

IFvif) = arg^gR{Fy(y) ^ f}, fe [0, 1], 
IFxiYif,y) = '<^rg,^AFxiY{^\y) = /}, fe [0,1]. (3) 

First, the y-direction position of the i*^ particle is di- 
rectly taken from the marginal ICDF as yf = IFY{{i — 
0.5) /N),i = 1,2,- •• where N is the number of 

particles and the superscript s indicates that yf is a 
sorted, monotonically ascending sequence. Then, corre- 
sponding to each y| generated in this manner, the x- 
direction position of the i*'* particle is determined as 
Xi = IFx\Y{Ui,yi), where Ui are i.i.d. real random 
variables with uniform distribution over [0,1]. We only 
have a finite number of conditional CDFs, and we will 
assume that they are smooth in y; therefore, for each 
particular working y| we employ the conditional CDFs 
available in its neighborhood (e.g. the closest one, or 
possibly an interpolation of the closest ones). Only a 
few conditional ICDFs are needed if the CDF is suffi- 
ciently smooth. For example, if M {M <C N) conditional 
ICDFs are needed, then these ICDFs, /Fx|y(/, y£), fc = 
1,2,- •• ,M,f S [0,1], can be chosen such that ?/£ ~ 



fc-i)-int(Af/Af)+int(Af/2Af) 



, where int(^),^ S R is the 



maximum integer not greater than ^. 

2.2. Fine-Scale Evolution. We work with mul- 
tiscale dynamical systems for which we have fine-scale 
evolution models available. For noninteracting particles, 
letting X{t) = {x{t),y{t))'^ £ denote the fine-scale 
model state, consisting of particle positions at time t, 
the discrete dynamics for X are given by 



Xk+i='f{Xk,Vk;At;X), k^O,l, 



(4) 



where ijk G denotes an external stochastic driving 
force and A is a set of (constant) parameters. In the 
class of problems envisioned here, the models \E'(-) involve 
microscale simulators of Brownian motion, kinetic Monte 
Carlo simulation or molecular dynamics. 

2.3. Restriction. Let the position of the i*'' parti- 
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particle positions in two directions be denoted as {Xj} 
and {yfn}, respectively. Then a two-dimensional mesh 
can be formed with each grid point having a coordinate 
{Xj,y!„),j, TO = 1, 2, • • • , iV. For each point (xj, ?/f„), the 
number, Nf, of particles whose x and y direction po- 
sitions satisfy Xi < Xj and yi < y!^^, respectively, is 
counted and the CDF at this grid point evaluated as 

PxYixj,ym) = "^^^ CDF of particle positions 

can thus be obtained. This is only one of several possible 
restriction methods; other restriction approaches can be 
found ic^^. 

Assuming the CDF FxY{x,y) to be differentiable, 
Eqn. (QJ or (O lead to the following fomula for the con- 
ditional CDF, 



dFx 



X\Y 



dy 



,y) 



dFY 
dy 



(y) 



(5) 



Continuing the example at the end of Section 2.1, the 
evaluation of the conditional ICDFs, (/, y^), fc = 

1, 2, • • • , M, / £ [0, 1] , requires the availability of the con- 
ditional CDFs, Fx\Y{x\y^)- These can be numerically 
approximated from equation jSJl as. 



Fx\Y{x',\yl) 



FxYjx^j.ypJ - FxY{x1,yl^) 
FxY{x%,y^J - FxY{x%,y;J' 



(6) 



where pi and p2 can be chosen as (fc — 1) • int(iV/A/) + 
1 and k ■ int(iV/M), respectively. Once the conditional 
CDF Fx\Y{x\y'j^) is available numerically, the conditional 
ICDF IFx\Y{f,yk) can be numerically evaluated as in 
the case of one-dimensional observables. 

3. Coarse Projective Integration (CPI) for Mul- 
tidimensional Random Particle Systems 

Coarse projective integrators (CPI) typically consist of 
four steps (Fig. [SJ. At first, coarse observables are iden- 
tified to which the fine-scale model states are restricted. 
The coarse observables used in our context consist of the 
marginal and (finitely many) conditional ICDFs of the 
microscale particle positions, as described in the previ- 
ous section; in particular, we use a finite number of ex- 
pansion coefficients (in some appropriate basis) of this 
marginal and these conditional ICDFs. The particle po- 
sitions can be generated through the lifting procedure de- 
scribed in Section 2.1 once an initial condition for these 
ICDFs have been specified. We represent the lifting op- 
erator by /i; a mapping from the coarse observables IF 
(ICDFs) to the microscopic descriptors X (particle posi- 
tions). The second and third steps are the fine-scale evo- 
lution and restriction mentioned in sections 2.2 and 2.3, 
respectively. The restriction operator M. is a, mapping 
from the microscopic descriptors X to the coarse observ- 
ables IF, i.e., IF = MX. Evidently, the operators M. 
and ^ satisfy the property Mfi = I (modulo roundoff 

prrnr^ Alnncr wif.h rpstrirf.inn rnmps f.hp psf.impii.inn nf 
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the coarse-scale time derivatives of the observables (the 
marginal and ICDF coefficients). 

The last step is the projection step in time - the tem- 
poral evolution of our representation of the coarse-scale 
observables. This step is templated on continuum nu- 
merical integration techniques - for coarse forward Euler 
it is simple linear extrapolation of the coarse observables 
in time, although more sophisticated and even implicit 
techniques can be (and have been) usedi^*2&. 

Let the coarse-scale observables at time t consist of 
M + 1 ICDFs, /F,,t,i = ,M + 1, of which the 

first one is the marginal ICDF, IFyif), and the remain- 
ing ones are the conditional CDFs, /-Fx|y(./, 2/fc), at y- 
dircction positions y'j.,k ~ 1, 2, • • • , A/. Let a basis of the 
coarse-scale subspace, which can be specified globallj*^ 
or localljii, be denoted by {9q}. Then the projections, 
Pi^q^t, of ICDFs onto the basis can be computed by, 

P^,q,t^iIF,,t,eq), 9 = 0,1,. ..,P t=l,2,..-,n, 

(7) 

where the inner product is in L2- 

The projective integration step over a coarse-scale time 
step T can be formally written as 

Pi,q.n+T = L{(3i^q,l^l, ■ ' ' :Pi,q,n), (8) 

where L(-) is an operator based on (templated on) tradi- 
tional continuum numerical integration schemes, (3i_q_n+T 
is the coefficient of the q*^ mode of the z*'* ICDF imme- 
diately after the temporal projection step and Pi^q^t,t = 
l + l, - ■ ■ , n is the coefficient of the q*'' mode correspond- 
ing to the i*'* fine-scale time step prior to the projective 
step. 

Immediately after the projective sep, the new ICDFs, 
IFi^n+T, based on the new coefficients I3i_q^n+T, are con- 
structed as, 

p 

IFi^n+T ^"^^Pi^q.-n+T^q ■ (9) 

q=Q 

New fine-scale model states can then be "lifted" from the 
ICDFs, IF,^n+T- 

Using the above steps, the procedure for applying CPI 
to multidimensional random particle systems can be sum- 
marized as follows: 

1. Generate fine-scale model state(s) consistent with 
the coarse-scale description given by the particle 
ICDFs. The marginal and conditional ICDFs can 
be obtained through equation |SJ) if an analytical 
two-dimensional CDF is given as the initial condi- 
tion. 

2. Let the fine-scale model state evolve according to 

f.Vip rlisrrpf.p Hvnamippi.l Tnnr](=>l 



3. Generate ICDFs at some successive fine-scale time 
steps. 

4. Project the ICDFs onto an appropriate basis (equa- 
tion O), and estimate the temporal derivatives of 
coefficients of the dominant modes. 

5. Extrapolate (project forward in time) coefficients of 
the dominant modes over a large coarse-scale time 
interval T (equation (jSJ), reconstruct the ICDFs 
(equation ©) and go back to step 1. 

Usually one wants to report the multidimensional CDFs, 
and they can be generated numerically along with the 
ICDFs using the approach in Section 2.3 (although only 
the ICDFs are taken as coarse-scale states involved in the 
CPI method). The mesh size for numerically comput- 
ing the CDF can be set larger than maxi<i<jv_i(a;f_|_]^ — 
a;f,y|_|_i — 2/|), in an attempt to alleviate fluctuation- 
related problems in the estimation step; variance reduc- 
tion schemes (multiple realizations of the simulation, or 
more sophisticated approaches) may become necessary 
for this purpose. 

4. Coarse Dynamic Renormalization (CDR) for 
Multidimensional Random Particle Systems 

For multiscale systems of practical interest, if the 
PDEs in the macroscopic level are scale invariant, they 
may possibly possess self-similar solutionsi^. The anal- 
ogy with traveling wave solutions for problems with 
translational invariancc is instructive: approaching a 
traveling wave in a co-traveling frame appears like the 
approach to a stationary state. Similarly, approach- 
ing a self-similar solution in a dynamically renormalized 
(co-exploding or co-collapsing) frame, appears like the 
approach to a stationary state. Dynamic renormaliza- 
tion procedures have been used to investigate self-similar 
systemsiSiiLiS; recently a template-based approach for 
studying the dynamics of problems with translational 
symmetrjiSS has been extended to study the dynam- 
ics of problems with scale invariancoiSiSiiSSiS^. When 
macroscopic scale- invariant PDEs are explicitly available, 
template conditions can be used to derive dynamical 
equations (termed "MN-dynamics" ) for the rescaled self- 
similar solutions and similarity exponents^!. The idea of 
employing template conditions can also be used to obtain 
renormalized self-similar macroscale solutions and simi- 
larity exponents for multiscale systems whose coarse-level 
PDEs are not explicitly knowniS. The number of tem- 
plate conditions depends on how many rescaling variables 
arc needed to renormalizc the physical solutions. 

Consider a PDE in the form of 



^ = D.y{F), (10) 
where F{x,y^t) is a CDF of particle positions which do 

Tint rnllprtivplv trp.nslp.tp in tVip .cinnr'A rlnmain TtViA rnsp nf 
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FIG. 2: The evolution of the coarse and fine-scale model states. 



joint scale and translational invariance can also be simply 
treated, seoS^). 

The differential operator D^^y is such that there exist 
constants p and a such that 



X 



= JL 

(11) 

for any real function /, real value A > and coordi- 
nate [x, y) (there is no amplitude rescaling since this is 
a CDF). If a self-similar solution F{x,y,t) exists, it can 
be written as 



Fix,y,t)^Ui 



(cs)" ' (cs)"P 



;c). 



(12) 



where c is a constant parametrizing the family of self- 
similar shapes, s = t — to, {to is the blowup time for 
problems with finite time singularities) and t > to- Sub- 
stituting H12(l into H1Q(I . we have 



-1, 



and U satisfies the PDE, 



-auUu - apvU„ = c ^Duv{U), 



(13) 



(14) 



where u = x/{cs)°', v = y/(cs)°'P. For the operator D^y 
that satisfies Eqn. idJ, the constant a is determined 
by Dj;y itself. Hence the similarity exponent a can be 
calculated by Ean. lfT^ . 

If the macroscopic equation H10() is not explicitly 
known, one cannot analytically obtain the exponents 
p and a; numerical computations are needed to deter- 
mine these constants - and thus to test the scale invari- 
ance of the operator - before locating the self-similar 
solutions themselves. For an operator D^y that satis- 
fies the (unknown) equation (jll|l . the constants p and 
a can be obtained using a black box simulator of the 
equation as follows: Since the unknown Eqn. (|ll|l is 
valid for any coordinate {x, y) and real function /, let 
/ be a test function (we choose it here for convenience 

t.n V»=> PYnnnAntip.l in snarA^ Inr.-, Hi \ = A AP^ fi_Tirl 



{x2,y2) = {U2A, V2AP), where A is arbitrarily chosen as a 
positive real value. This would imply that (choosing two 
points in space) the following two relations hold: 



D.v{f{^. ^)){xuyi) = A'^D^,,{f{u, v)){u^,vi), 

D.yifij, J^))i^2, y2) - A'^D^,{f{u, V)){U2,V2), 

where {ui,vi) and {u2,V2) are two distinct coordinates. 
Comparing the above two equations, we have 



Dxy{f{% -J^)){X2,V2) Duvifiu, v)){u2, V2)' 



(15) 



Therefore, the constant p is the solution to Eqn. H15|l . 
We rewrite Ean. ljlSII as 



. X y 



Dxy{f{j,^)){x2,y2) 

D.y{f{%i,)){x,,y^) 
Duv{fiu,v)){ui,vi) 



(16) 



Duv{f{u,v)){u2,V2) = 0. 



The constant p can then be solved for, employing an 
in principle arbitrary test function /, and using New- 
ton's method. Since we assumed that the operator D^y 
is not explicitly available, Dxy{f{x,y)) can be estimated 
by running the micro-simulator for short time "bursts" 
and numerically obtaining the derivative The con- 
stant a is calculated by 



log. 



D.y{f{%i,)){x,,y,) 
Duv{f{u,v)){ui,vi) 



(17) 



once p is obtained. Clearly, other test functions and con- 
ditions evaluated at other points can be used; care must 
be taken also to ensure the finitcness of the estimated 
quantities. 

Given p, to determine the self-similar shape of the so- 
lution, we consider the general scaling 



F{x,y,t) = 



V 



(18) 



'A{t)' A{t)p 

where A{t) is an unknown function. The PDE becomes 
pAt 



At 

LUt rULUu 



A'^Duvi^). 



(19) 
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Evidently, U and lu are both renormalized CDFs. 
By comparing equations 114|l and (|19|l . we have 



(20) 



Also by comparing equations l|12(l and p8(l . we have 
limt^oo -4 



limt^oo (cs)*^ 



(21) 



The above equations and 1(211 together with itT^ 
lead to 



it-to)At 
a = hm ; . 

t^oo A 



(22) 



Therefore, the value for a can be calculated once values 
of A(t) are obtained in the long-time limit (i.e., after uj 
reaches the steady state). Indeed, let ti and t2 be distinct 
times after a; reaches the steady state, then by (|22|) . 



t2-ti 



Ait2) _ Ajti) 
At{t2) At{ti) 



(23) 



A single template condition is required to solve for both 
u)(u,v,t) and A{t) at every time step. In this paper, the 
template is chosen to be 

w(e,oo,t)=m, e < 0, < m < 0.5, (24) 

where e and m are both constants. The template con- 
dition has the following physical meaning: the rescaled 
marginal CDF lou always has the same value m at the 
u-direction coordinate e for all time t. Applying this 
template to Ean. (|19|l and assuming ^(e,f,t) decays ex- 
ponentially as u — > cx), we have 

Atef^(e,oo,i)+A'^+ii?™(w)(e,oo,t) = 0. (25) 
ou 

Equations (|19(l and H25(l can be coupled to solve for the 
rescaled CDF w and rescaling variable A if the operator 
Dxy is explicitly known. As the time t — > oo, a; may ap- 
proach a steady state, which is then a stable self-similar 
shape for the solutions to Ean. llO|) . In general cases, the 
macroscale equation for the CDF of the particle positions 
may not be explicitly available. However, the template- 
based approach can still be used to renormalize the CDF 
evolved via microscale models and rescaling variables are 
obtained during the course of renormalization. In these 
cases, to express the CDFs, the marginal and conditional 
ICDFs are used again as macroscopic observables, and 
their (discretized) projections over an orthonormal basis 
are again used to numerically characterize them through 
a finite number of coefficients. Based on the coarse time- 
stepper, the procedure for the coarse renormalization is 
schematically depicted in Fig. |2|and consists of the fol- 

Inwincr sf.Ans 



1. Generate the marginal and conditional ICDFs 
IFi^t ,z = l,---,M-|-l according to the initial CDF 
using equation © or © or according to the coeffi- 
cients, of dominant modes of the ICDFs using 
equation 

2. Generate particle positions using the ICDFs using 
the lifting procedure in the coarse time-stepper. 

3. Evolve particle positions over a time interval T' us- 
ing a fine-scale model (0J. 

4. Obtain ICDFs from particle positions using the re- 
striction procedure in the coarse time-stepper. 

5. Rescale the marginal ICDF according to the tem- 
plate condition and obtain the rescaling variable A. 
We then rescale the conditional ICDFs by a fac- 
tor of AP. This step can be justified by Ean. (|18|l . 
Indeed, obtaining the rescaled solution ujk+i from 
ujk via the dynamics (|19|1 and (|25|l is equivalent to 
starting from the initial condition uJk, running the 
original dynamics for a while to get i^fc+i, factoring 
out the rescaling variable A, and rescaling F^+i in 
scales of A and A^ respectively in x and y direc- 
tions. 

6. Project the rescaled ICDFs onto the orthonormal 
basis and obtain the coefficients of leading modes 
using equation {T)). Go back to step 2. 

The above procedure can be viewed as an iterative al- 
gorithm to solve the fixed point of a nonlinear operator 
^T' , written as. 



/3 = <I>T'(/3). 



(26) 



This fixed point can be written in component form as, 
Pl^J^l,--- ,M+l,q = 0,--- ,P^OY 

/3 = (/^l.O' /^l.l I ■ ■ ■ : Pl,P: ^2,0: ^2,1: ' ' ' 1/^2, Pi 

where the superscript r refers to the fact that this is 
the renormalized self-similar shape. These coefficients 
correspond to the renormalized self-similar ICDFs and 
CDF of the multidimensional particle system. Equation 
(|26|l may be solved using any numerical algorithm such 
as direct iteration or matrix-free (Krylov-subspace based) 
implementations of Newton's method^S. 

5. Self-Similar and Asymptotically Self-Similar 
Dynamics of Brov^rnian Particles in a Couette 
Flow 

We will use CPI and CDR algorithms to study two- 
dimensional Brownian models of particle dispersion in a 
Couette flow^. In this section, a particle system with 
self-similar dynamic evolution and a system with asymp- 

f.nf.irn.llv S(=>lf-siTnilp»r pvolntinn pta AvnlnrpH rpsnpptivAK;- 
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5.1. A Particle System with Self-Similar Dy- 
namics. Let X{t) and Y{t) represent particle positions 
in X and y directions respectively at time t on an infinite 
two-dimensional spatial domain. The particle positions 
in the two directions evolve in this model governed by 
the following dynamics: 



dX{t) = DdWxit), dY{t) = Xdt, 



(27) 



where Wx (t) is a Wiener processes^ and D is the diffu- 
sion coefficient. The discretized dynamics of (|27|l is given 



Xk+i =Xk + DTjx,kVAt, Yk+i =Yk+ XkAt, (28) 

where rix,k ai'c i.i.d. standard Gaussian random vari- 
ables. 

The dynamics (|27|) represent the motion of particles 
which only diffuse in the .x-direction in a Couette flow. 
It can be shown that the coarse-scale dynamics for the 
PDF. PxY{x,y,t), of a particle position, corresponding 
to the fine-scale dynamics (|27|) . is governed by the fol- 
lowing equations^ 



dPxY , dPxY D^d^PxY 



dt 



dy 



(29) 



where Pxy is assumed to be 2nd-order differentiable. 
Hence the dynamics for the CDF, FxY{x,y,t), associ- 
ated with (|29|l is given by 



dF 



XY 



dF 



XY 



dF- 



XY 



dx 



d^FxY 



(30) 



dt ' ~ dy dy " 2 dx'^ 

In the above equation, the operator D^y is written as 



d 



d^ 



^^^ = -^d-y^l^^d-y'^'--dx^- ^''^ 

This operator satisfies the scale invariance property 
for constant values p = 3 and a = — 2. The ana- 

K/-tippiI SAlf-similpiT snlntinn tct R,nn H'/Qh is insnirpd Vi\^^ 



(See Appendix A) 



PxY{x,y,t) 



7:D^t-toy 



6(y-0.53;(t-to))"' 



^ 2D''' (t- 



(32) 

where tg is the blowup time (backward in time), which 
then gives the self-similar solution to H3Q|1 . 



FxY{x,y,t) = 



73 



, 6tH-0.5x(t-to))^ ^2 

^ i^^(t-'o)^ ^^'^Ht-to)>dydx. (33) 



OO — oo 



Let u' 



and v' 



(c(t-to))V 

FxY{x,y,t) = Fuv{u',v') 



(c(t-to))3/2' 

V3 



then 



7rL»2/c2 



e '^^^-'dvdu. (34) 



OO — oo 



Hence for the integro-differential equation H3U|1 . the 
similarity exponent in (|12|l is a = 1/2. For the CDF 
in H34|) . its standard deviations (std.'s) in two directions 
and correlation are ax = D/c^^^, (Jy = D/iV^c^^"^), and 
Pxy = V3/2, respectively. 

5.2. A Particle System with Asymptotically 
Self-Similar Dynamics. We now consider Brownian 
particles in a Couette flow that diffuse in both spatial 
directions. The microscopic-level evolution equation for 
the particle positions is given by 

dX(t) ^ DdWxit), dY{t) = Xdt + DdWY{t), (35) 

where Wx{t) and WY{t) are independent Wiener pro- 
cesses. The discretized dynamics of p5|l are given by^^ 



Xk+i = Xk + Dr]x,ky/Ai, 
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where r]x,k and r]Y.k arc i.i.d. standard Gaussian random 
variables. 

The coarsc-scalc PDF of the particle positions is gov- 
erned bji2^ 



dPxY , dPxY D^d'^PxY , D^d^PxY 



dt 



dy 2 dx'^ 2 dy'^ 



(37) 



Accordingly, the dynamics for the CDF, Fxy{x, y, t), are 
given by 



dFxY , dFxY 
- X 



dF 



XY 



dx 



dt dy dy 

d^FxY d^FxY 
~2 dx^ ^ ~2 dy"^ ' 

In Eqn. H38|l . the operator D^y is written as 



(38) 



d r d , d'' ^2 

The above operator does not exactly satisfy the scale 
invariance property Hll|l for any function /, positive real 
value A and coordinate {x,y). However, as the scale of 
the function / or the value of A become sufficiently large, 
there may exist p and a such that is approximately 
satisfied. At those limits for / and/or A, values of p and a 
that approximately satisfy approach the values that 
exactly satisfy the scale invariance property in the true 
self-similar case. 

The asymptotically self-similar solution to the equa- 
tion is given bji^i 



PxY{x,y,t) 



1 



2-KD^{t-t^){l + [t-t^Y/12Yl^ 



(B-0.5i(t-to))^ 



2-D2(t-t(,){l + (t-to)2/12) 2D^(t-tQ) 



(40) 



which then provides the asymptotically self-similar solu- 
tion to 

1 

FxY{x,y,t) = 



2TTD^it-to)ii + {t-to)yi2y/^ 

(y-0.5x(t-to))2 



-D2(t-to)(i+{t-f 0)^/12) '^o^(t-to)' dydx. 



— oo — oo 



(41) 



Let u' 



TTT2 and v' 



then 



FxY{x,y,t) = Fuv{u ,v' ,t) = 
Vi{t - to) 

TTD^{{t-tof + 12)^/^1 

r' . 6(^-0.5u/c)2(t-to)2 „2 

/ e ^ c2((t-to)2+i2)/c:i ^^^-'dvdu. (42) 

-oo — oo 

In the long-time limit, Fijy{u' , v' , t) has a steady-state 
form, which is the same as that given by (|34|l in the self- 
similar rpisp 



6. Numerical Examples 

In what follows, direct particle simulations are im- 
plemented to accelerate the numerical evolution of the 
CDFs via CPI, and to locate self-similar CDFs via CDR 
for the particle systems in Section 5. The approximate 
(asymptotic) scale invariance of the macroscale differen- 
tial operator D^y for the asymptotically self-similar par- 
ticle system is also examined. The fixed-point algorithm 
in Section 4 is utilized to solve for the long-time steady- 
state shape of the CDF for the asymptotically self-similar 
system. The diffusion coefficient D and simulation time 
step Ai in the fine-scale model are set to S.Ocm/s^/^ and 
0.01s, respectively. An ensemble of 2000 [N = 2000) par- 
ticles is used in the fine-scale simulations except where 
otherwise indicated. 

Simulation 1: Direct Simulation of the Self- 
Similar Particle System 

In this simulation, the initial fine-scale particle posi- 
tions are chosen to follow a uniform distribution over 
the square domain (—10cm, 10cm) x (—10cm, 10cm). An 
ensemble of 2000 particles whose distribution is consis- 
tent with the coarse-scale initial conditions are evolved 
directly using H28|) and used to construct true evolved 
coarse-scale CDFs. Particle positions are recorded at 
the time 300At, 600At and 900At, respectively. Two- 
dimensional CDFs are numerically computed using the 
procedure in Section 2.3 and plotted in Fig. 01 The num- 
ber of grid points used to compute the CDFs is 1681. 

Simulation 2: Coarse Projective Integration of 
the Self-Similar Particle System 

We now use CPI to accelerate the evolution of the 
coarse-scale observables. The coarse-scale initial condi- 
tion is the same as that in Simulation 1. Particles are 
evolved for an initial block of 10 {I — 10) fine-scale time 
steps and then again for another block of 10 {n — I = 10) 
fine-scale time steps. At each of the latter 10 steps, the 
marginal ICDF and 20 (M = 20) conditional ICDFs of 
particle positions are formed; the time series of the coef- 
ficients of their leading modes is linearly extrapolated 
(with a slope estimated through least-squares) over a 
time interval equal to 10 (T = 10) fine-scale time steps. 
The basis onto which the ICDFs are projected consists 
here of shifted Legendre polynomials of order up to and 
including 5 (P = 5). Since the ICDFs are anti-symmetric 
with respect to the axis / = 0.5, the coefficients of the 
2nd and 4th modes vanish. Hence, only 84 coefficients 
need to be extrapolated. At the end of the extrapola- 
tion, the ICDFs are reconstructed and particle positions, 
generated according to these ICDFs, are simulated again. 
At time 300At, 600Ai and 900At, particle positions are 
recorded and CDFs plotted in Fig. [S] The number of 
grid points used to compute CDFs is again 1681. The 
cross sections, Fxy{s, s, t), of CDFs in Simulation 1 and 
2 are compared in Fig. El which shows an excellent visual 
match between the true CDFs and those computed from 

f.ViA r^PT p.lcrnrif.tim 
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No. of Iterations 


P 


a 





5.0 


-3.34959 


1 


2.80093 


-1.88311 


2 


2.99246 


-2.03967 


3 


3.00106 


-2.04212 


4 


3.00370 


-2.05468 


5 


2.99592 


-2.03927 


6 


2.99659 


-2.03914 


7 


2.99753 


-2.04017 


8 


2.99831 


-2.04189 


TABLE I: Iterative evaluation of the constants p and a. 


No. of Iterations 


P 


a 





5.0 


-3.38239 


1 


3.19621 


-2.25167 


2 


2.97334 


-2.06146 


3 


2.99987 


-2.08557 


4 


2.99574 


-2.08444 


5 


2.99347 


-2.08425 


6 


2.99747 


-2.08345 


7 


2.99247 


-2.08073 


8 


2.99833 


-2.08601 



TABLE II: Iterative values for the constants p and a for a 
different set of algorithm parameters (see text). 



Simulation 3: Coarse Dynamic Renormalization 
of the Self-Similar Particle System 

Pretending that the macroscale equation H3()(l is not 
available, we now directly use the microscale simulator 
(|28|l to compute the constants p and a, the similarity 
exponent a and the macroscopic self-similar solution. 

In the approach provided in Section 4, Newton's 
method is used to solve the equation H17|) for p. The test 
function / was first chosen as a 2-dimensional joint Gaus- 
sian distribution function, f{x,y) = l/16N{x/4)N{y/4), 
where N{x) and N{y) arc standard Gaussian distribu- 
tions. We select the value of the positive real number 
A = 2.0. The two coordinates and (u2,V2) in 

(|17|l arc chosen as (—2,-2) and (3,3), respectively. To 
reduce fluctuations of values for the operator D^y, 9000 
particles are used and 500 replica copies of values for D^y 
and Duv are averaged in the computation. 

Starting from the initial value po = 5.0, iterative values 
for p are stabilized at 3.0 after 3 iterations. Accordingly, 
the converged value for a is —2.0 (Table P). 

For verification, we choose another set of parameters 
to determine p and a: f{x,y) = l/25N{x/5)N{y/5), 
A = 2.5, {ui,vi) = (-3,-3), and (^2,1^2) = (4,4). 
Again, starting from the initial value po = 5.0, iterative 
evaluations of p are stabilize at 3.0 after 3 iterations and 
the converged value for a is again around -2.0 (TableHHl. 

We can therefore conclude that there exist p = 3.0 and 
a = —2.0 such that the unavailable differential operator 
Dxy corresponding to the microsimulator H27|) possesses 
the scale invariance property (|ll|l . Accordingly, the sim- 

ilnritv A-ynnriATit is n/ = ^ lr\M E,nn 111 'Ah 



t{sec) A{t) At{t) 

~ 1.00000 

1 1.10268 0.10268 

3 1.27793 0.08763 



TABLE III; The rescaling variable A{t). 



The template condition for the x direction is chosen 
to be a;(— 2.832, oo, t) = 0.4, i.e., the u-coordinate corre- 
sponding to the renormalized marginal CDF uju = 0.4 
always has the same value, —2.832cm. The constant c in 
the analytical solution H34|) is obtained as c = 0.2sec^^ 
based on our template. The corresponding std.'s for 
the analytical self-similar shape are ax = 5\/5cm and 
try = 25v'T5/3cm, respectively. 

The CDF corresponding to a uniform distribu- 
tion of particle positions over the space domain 
(—10cm, 10cm) x (— lOcm, lOcm) is used as the initial 
condition. Direct iteration is used to solve for the fixed 
point of equation H26|) . The time interval T' is 100 At. 
The number of conditional ICDFs is 20 (M = 20) and 
the basis for the ICDFs is again the shifted Legendre 
polynomials of order up to and including 5 (P = 5). 
In this simulation, 100 copies of ensemble particle po- 
sitions are generated according to the mode coefficients 
of the ICDFs at the beginning of each iteration and let 
to evolve. The mode coefficients at the end of each it- 
eration are obtained by averaging over these 100 replica 
copies. After the 2nd, 4th and 6th iterations, renormal- 
ized mode coefficients of the ICDFs are used to generate 
particle positions, out of which the CDFs are computed 
and plotted respectively in Fig. Fig. ISjalso compares 
the cross sections u!{u,u,t) of true CDFs and renormal- 
ized CDFs in this simulation. Clearly, the renormalized 
solutions quickly approach the self-similar steady state. 

To validate the computation of the self-similar solution 
shape, the std.'s and correlation of the computed shape 
arc compared with those of the known analytical solution. 
The std.'s and correlations of the rescaled CDFs are cal- 
culated via the ensemble particle positions corresponding 
to these CDFs. The comparison is shown in Fig. |51 where 
curves in Case 1 represent results obtained using this 
template condition and time interval. The std.'s and cor- 
relations of the rescaled CDFs approach those of the an- 
alytical self-similar shape, which means that the rescaled 
CDF coincides eventually with a member in the family 
of theoretical self-similar shapes expressed by Equation 

As the renormalized CDF cu reaches its steady state, 
we can set this CDF as the initial condition and evolve 
the microscale dynamics (|28|) for two more loops with 
h = lOOAt and <2 = 300At. The rescaling variable A{t) 
is listed in Table |llll Note that A{t) = 1 at t = 0. By 
Equation H23|l . the similarity exponent a is approximated 
as 0.520, within 4% of the theoretical value 1/2. 

In the following, the effect of variation of templates 

anH Avnlntinn timps nf thp. fiYPrl--nnint nnAratnr <I>^/ nn 
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No. of Iterations 


P 


a 





5.0 


-3.87764 


1 


4.14589 


-3.32940 


2 


3.81151 


-3.09345 


3 


3.78553 


-3.08460 


4 


3.79408 


-3.08433 


5 


3.78576 


-3.08486 


6 


3.79758 


-3.08472 


7 


3.78701 


-3.09319 


8 


3.79517 


-3.08471 



TABLE IV: Iteratively computed values of p and a for Pa- 
rameter set 1 (see text). 

the computed renormalized self-similar shapes will be ex- 
amined. For self-similar systems, we can see from Eqn. 
()19|1 that, as the system reaches the steady state, the 
rescalcd shape of CDFs will remain the same, irrespec- 
tive of changes in the evolution time T' . Also, the steady- 
state CDF shapes will coincide with members in the fam- 
ily of self-similar solutions prescribed by Eqn. H14|l no 
matter what the scale of the template is. 

We choose four cases of the template condition and 
evolution time including the one above: 

1. cj(-2.832,oo,t) = 0.4, T' = lOOAt; 

2. cj(-2.832,oo,t) = 0.4, T' = 200At; 

3. cj(-0.283,oo,t) = 0.4, T' = lOOAi; 

4. cj(-0.283,oo,t) = 0.4, T' = 200Ai. 

The iterative values of std.'s and correlation for the 
four cases are shown in figures El and 1101 Comparison 
with theoretical calculations shows that variation of tem- 
plates and evolution times indeed does not cause devia- 
tion of the converged rescaled CDF from the family of 
self-similar solutions. 

Simulation 4: Coarse Dynamic Renormalization 
of the Asymptotically Self-Similar Particle Dy- 
namics 

For the particle system in Section 5.2, the procedure 
in Section 4 is used to check if its macroscopic differen- 
tial operator D^y possesses the scale invariancc property 
(|ll|l . The two parameter sets in Simulation 3 are used 
here. 

1. Parameter set 1: f{x,y) = l/l%N{x/A)N{y/4), 
A = 2.0, K, fi) = (-2, -2), (U2, V2) - (3, 3); 

2. Parameter set 2: f{x,y) = l/25iV(a;/5)7V(?;/5), 
A - 2.5, = (-3,-3), {U2,V2) = (4,4). 

Newton's method is utilized again to solve for p and a 
for each parameter set. The iteratively computed values 
of p and a are listed in tables IIVI and El 

It can be seen that now converged values of p and 

n. dn imrii with t.ViA tprrmlntp spa.lp anrl valnp nf A As 



No. of Iterations 


P 


a 





5.0 


-4.23757 


1 


4.47814 


-3.85506 


2 


2.58903 


-1.95833 


3 


4.01019 


-3.38693 


4 


3.64906 


-2.96049 


5 


3.35744 


-2.69786 


6 


3.39829 


-2.74037 


7 


3.40776 


-2.75770 


8 


3.41325 


-2.74170 


9 


3.41298 


-2.75793 



TABLE V: Iteratively computed values of p and a for Param- 
eter set 2 (see text). 



the template scale and A increase, values of p and a ap- 
proach those in the self-similar case. Pretending that the 
macroscale equation is not explicitly known, we may sus- 
pect that the particle system exhibits asymptotically self- 
similar dynamics. Therefore, for "asymptotically large 
enough" template conditions, the operator D^y still ap- 
proximately possesses the scale invariance property (|ll|l ; 
as we did in the self-similar case, we can use a fixed 
point algorithm to find a long-time steady state for the 
asymptotically self-similar solution. In analogy to the 
self-similar case, the evolution time interval in the fixed- 
point operator does not affect the converged shape for 
such large enough scales. 

The four template condition and evolution time cases 
in Simulation 3 are used in the fixed point algorithm to 
verify the above assertions. The value p is set to 3.0, the 
same value as that in Simulation 3. The microsimulator 
is the discretized dynamics (|36|) . Iterative values of std.'s 
and correlation for the four cases are shown in figures [TTI 
andll2l As can be seen, for large template conditions, the 
length of evolution time T' does not affect the converged 
values of std.'s and correlation for distribution of particle 
positions. Yet for small templates, the effect of T' on 
the y-direction std. and correlation is evident. As T' 
increases from lOOA^ to 200A<, the converged values of 
?/-direction std. and correlation approach their expected 
theoretical values. 

7. Conclusions and Remarks 

We presented an equation-free computational ap- 
proach, based on using marginal and conditional ICDFs 
as coarse-scale observables, for the computer-assisted 
study of multidimensional random particle system dy- 
namics. Coarse projective integration employing this 
time-stepper can be applied to accelerate the computa- 
tional evolution of particle CDF computations; the ap- 
proach targets multidimensional particle systems whose 
coarse-scale models are not explicitly available. Coarse 
dynamic renormalization can also be used to analyze 
particle systems with self-similar or asymptotically self- 
similar coarse-grained evolution dynamics, and to obtain 
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tions. 

The examples in this paper are admittedly simple, yet 
they illustrate the computational approaches in a con- 
text where the results can be validated; we hope that 
the type of multiscale algorithms presented here may be 
useful in more complicated situations (e.g., particles mix- 
ing in time-dependent velocity fields) if the macroscop- 
ical dynamic arc effectively self-similar. Another pos- 
sible application of such equation-free approaches is in 
cases where even the coarse-scale observables are charac- 
terized by uncertainty/stochasticity. Polynomial chaos 
observables have been used in the solution of explicit 
macroscale PDEs for passive scalar transport, where the 
uncertainty enters through random initial conditions or 
boundary conditions^Si^. Such polynomial chaos observ- 
ables may be combined with the coarse-graining tech- 
niques presented here when no explicit coarse-scale de- 
scriptions of the particle system dynamics are available. 

Appendix A 

For Equation 129|) with the initial condition 

PxY{x,y,to) ^ S{x)S{y), 

(where to is the blowup time) it is shown in the following 
that the solution is 

■nD'^[t - toY 
We can see that 

dt 

, -2V3 V3 {x-yf 

18{y-0.5x{t^to))^ 6iy-0.5x{t-tn))x 

dy 

-V3x 12(7/- 0.5a-(t- to)) 

D_()_FxY_ ^ n2(t-to)S +2nHt-to)) 

2 9x2 2 

-4^3 x/3 {6{y) - 4x{t - to))^ ' 

T:D^{t-taY T:D^{t-toY' D^{t-t^Y 

Hence, 




FIG. 4: True CDFs at a sequence of time steps; top left: 
t = Q, top right: t = SOOAt, bottom left: t = 600Af, bottom 
right: t = 900At. 



dPxY ^ JPxY d^PxY 
dt dy 2 dx'^ 
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FIG. 5: CDFs at different time instances computed by the 
CPI method; top left: t = 0, top right: t = SOOAi, bottom 
left: t = 600At, bottom right: t = 900A£. 



FIG. 7: Renormalized CDFs in Simulation 3; top left: initial 
CDF, top right: renormalized CDF after 2nd iteration, bot- 
tom left: renormalized CDF after 4th iteration, bottom right: 
renormalized CDF after 6th iteration. 



true CDF, t-0 
true CDF, t-300dt 
true CDF, t-600dt 
true CDF, t-900dt 

CPI CDF, t-0 
CPI CDF, t-300dt 
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CPI CDF, t-900dt 
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true CDF, t-300dt 
true CDF, t-600dt 
true CDF, t=900dt 
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FIG. 6: Comparison between cross sections of true CDFs and 
CDFs computed by CPI. 



FIG. 8: Comparison between cross sections of true CDFs and 
rescaled CDFs in Simulation 3. 
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FIG. 9: Comparison for standard deviations and correlations 
of self-similar shapes in Simulation 3: Cases 1,2. 



FIG. 11: Comparison for standard deviations and correlations 
of converged shapes in Simulation 4: Cases 1,2. 
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FIG. 10: Comparison for standard deviations and correlations 
of self-similar shapes in Simulation 3: Cases 3,4. 



FIG. 12: Comparison for standard deviations and correlations 
of converged shapes in Simulation 4: Cases 3,4. 
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